Ειδικά Μαθήματα Βοτανικής
1 Βασικές προϋποθέσεις
Εννοείται ότι πριν ξεκινήσουμε να κάνουμε το οτιδήποτε, έχουμε δημιουργήσει ένα νέο project στο R-Studio, το όνομα του οποίου - για λόγους τους οποίους θα καταλάβετε αργότερα - ΔΕΝ πρέπει να περιέχει:
1. ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)
Το ίδιο ισχύει ΚΑΙ για το file path του εν λόγω φακέλου (π.χ., όχι ‘E:/Ειδικά Μαθήματα Βοτανικής/SDM,’ αλλά ’E:/Eidika_Mathimata_Botanikis/SDM).
Ένα άλλο κρίσιμο σημείο είναι το εξής: θα χρειαστεί να εγκαταστήσετε τα πακέτα που πρόκειται να χρησιμοποιήσετε σε αυτό το tutorial1.
2 Εγκατάσταση και φόρτωση βιβλιοθηκών
Όπως κάθε φορά, θα χρειαστεί να εγκαταστήσουμε και να φορτώσουμε τις απαραίτητες βιβλιοθήκες.
## ===========================================================================##
## Install the main packages
## ===========================================================================##
install.packages(c("tidyverse", "rgbif", "raster", "data.table", "rgeos", "rgdal",
"gridExtra", "dismo", "biomod2", "sf", "usdm", "speciesgeocodeR", "pacman", "spThin",
"CoordinateCleaner", "ggspatial", "scales", "scrubr", "mapr", "rasterVis", "ade4",
"viridis", "biogeo"), dependencies = T)
## ===========================================================================##Ας τις φορτώσουμε:
## ===========================================================================##
## Load them
## ===========================================================================##
pacman::p_load(tidyverse, rgbif, raster, usdm, rgeos, rgdal, usdm, biomod2, gridExtra,
countrycode, devtools, sf, spThin, biogeo, speciesgeocodeR, pacman, dismo, CoordinateCleaner,
ggspatial, scales, scrubr, mapr, rasterVis, data.table, ade4, viridis)
## ===========================================================================##3 Διανυσματικά δεδομένα
3.1 Σύνορα χωρών
Πρώτα θα κατεβάσουμε τις γεωγραφικές μεταβλητές για την χώρα που μας ενδιαφέρει (η Ελλάδα στην προκειμένη περίπτωση). Θα πρέπει να εισάγουμε (στην εντολή country στον κώδικα που παρατίθεται πιο κάτω) τον 3-ψήφιο ISO κωδικό για την χώρα μας που τον βρίσκουμε από αυτόν τον ιστότοπο ή με την εντολή ISO_countries <- getData("ISO3") %>% as.data.frame2. Εάν θέλουμε να περιορίσουμε την περιοχή μελέτης σε μια μικρότερη περιοχή, τότε θα πρέπει να χρησιμοποιήσουμε ένα άλλο λογισμικό, το Qgis, ώστε να “κόψουμε” την περιοχή που μας ενδιαφέρει. Προς το παρόν, δεν θα χρειαστεί να το κάνουμε αυτό, αλλά θα περιοριστούμε να κατεβάσουμε τα πλήρη δεδομένα για την Ελλάδα: .
## ===========================================================================##
## Download the data for Greece
## ===========================================================================##
## There are three levels available, you are free to explore them
Greece <- getData("GADM", country = "GRC", level = 3)
## ===========================================================================##Μπορούμε να αναπαραστήσουμε γραφικά τα δεδομένα αυτά.
## ===========================================================================##
## Plot Greece
## ===========================================================================##
plot(Greece)Είναι αρκετά εύκολο στην προκειμένη περίπτωση να περιορίσουμε την έκταση του πολυγώνου στην Κρήτη, με τις ακόλουθες εντολές:
## ===========================================================================##
## First we need to see what is inside the object we created
## ===========================================================================##
Greece## class : SpatialPolygonsDataFrame
## features : 326
## extent : 19.37236, 29.6457, 34.80069, 41.74801 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 16
## names : GID_0, NAME_0, GID_1, NAME_1, NL_NAME_1, GID_2, NAME_2, NL_NAME_2, GID_3, NAME_3, VARNAME_3, NL_NAME_3, TYPE_3, ENGTYPE_3, CC_3, ...
## min values : GRC, Greece, GRC.1_1, Aegean, <U+0386><U+03B8><U+03C9><U+03C2>, GRC.1.1_1, Athos, <U+0386><U+03B8><U+03C9><U+03C2>, GRC.1.1.1_1, Abdera, Abelokipi-Menemeni, <U+0386><U+03B8><U+03C9><U+03C2>, Dímos, Municipality, NA, ...
## max values : GRC, Greece, GRC.8_1, Thessaly and Central Greece, Tessa<U+03BB><U+03AF>a <U+03BA>a<U+03B9> Ste<U+03C1>e<U+03AC> <U+0395><U+03BB><U+03BB><U+03AC>da, GRC.8.2_1, West Macedonia, Tessa<U+03BB><U+03AF>a, GRC.8.2.9_1, Zografou, Zografos, Tessa<U+03BB><U+03BF><U+03BD><U+03AF><U+03BA><U+03B7><U+03C2>, Dímos, Municipality, NA, ...
names(Greece)## [1] "GID_0" "NAME_0" "GID_1" "NAME_1" "NL_NAME_1" "GID_2"
## [7] "NAME_2" "NL_NAME_2" "GID_3" "NAME_3" "VARNAME_3" "NL_NAME_3"
## [13] "TYPE_3" "ENGTYPE_3" "CC_3" "HASC_3"
Crete <- subset(Greece, NAME_1 == "Crete")## ===========================================================================##
## Plot Crete
## ===========================================================================##
plot(Crete)3.2 GBIF
Στο σημείο αυτό, είμαστε πλέον στην ευχάριστη θέση να μπορούμε να χρησιμοποιήσουμε τα δεδομένα θέσης για τα είδη τα οποία μας ενδιαφέρουν. Υπάρχουν δύο περιπτώσεις, όσον αφορά την απόκτηση τέτοιων δεδομένων:
- να τα έχουμε αποκτήσει πρωτογενώς (δηλαδή να έχουμε εργαστεί στο πεδίο και να έχουμε εντοπίσεις τις θέσεις εμφάνισης του εκάστοτε είδους και των υποπληθυσμών του)
- να τα έχουμε αποκτήσει δευτερογενώς (δηλαδή είτε να έχουμε χρησιμοποιήσει τις διαδικτυακές βάσεις ψηφιοποιημένων - και ενίοτε γεωαναφερμένων - δεδομένων GBIF και BIEN, είτε να έχουμε ψηφιοποιήσει χάρτες εμφάνισης θέσης3)
Το ιδανικό βεβαίως είναι να έχουμε αποκτήσει πρωτογενώς τα στοιχεία αυτά, καθότι στην περίπτωση αυτή, μπορούμε να είμαστε σίγουροι για τις γεωγραφικές συντεταγμένες των υποπληθυσμών του προς μελέτη είδους. Στην πλειονότητα των περιπτώσεων όμως, δεν έχουμε πρόσβαση σε τέτοια δεδομένα, οπότε αναγκαστικά καταφεύγουμε στην δεύτερη λύση και δη στις διαδικτυακές βάσεις δεδομένων. Ευτυχώς για εμάς, υπάρχουν κάποια πακέτα4 στην R, όπως το rgbif και το rbien, τα οποία αυτοματοποιούν την διαδικασία μεταφόρτωσης δεδομένων από τις διαδικτυακές αυτές βάσεις5.
Εμείς θα εργαστούμε με το πακέτο rgbif προς το παρόν.
Σήμερα – όπως και στην προηγούμενη εργαστηριακή άσκηση – θα ασχοληθούμε ένα εκ των έξι μονοτυπικών γενών της ελληνικής χλωρίδας, την Petromarula pinnata.
Fig 1. The endemic species Petromarula pinnata from Crete.
Πρώτα χρειάζεται να βεβαιωθούμε ότι το προς μελέτη είδος, υπάρχει όντως εντός της διαδικτυακής βάσης την οποία επιθυμούμε να χρησιμοποιήσουμε, προκειμένου να εξάγουμε δεδομένα. Η εντολή name_suggest() ερευνά όλα τα είδη (taxa καλύτερα) τα οποία μοιάζουν με το είδος το οποίο ψάχνουμε και κρατά μόνο το όνομα των taxa το οποίο ταιριάζει ακριβώς με αυτό το οποίο μας ενδιαφέρει6.
spp_petromarula <- name_suggest(q = "Petromarula pinnata", rank = "species", limit = 10000)
spp_petromarula## Records returned [2]
## No. unique hierarchies [0]
## Args [q=Petromarula pinnata, limit=10000, rank=species, fields1=key,
## fields2=canonicalName, fields3=rank]
## # A tibble: 2 x 3
## key canonicalName rank
## <int> <chr> <chr>
## 1 3167989 Petromarula pinnata SPECIES
## 2 11108006 Petromarula pinnata SPECIES
Στο βήμα αυτό αναγνωρίσαμε τον μοναδικό αριθμό αναφοράς (key) του προς μελέτη είδους και μόνο χρησιμοποιώντας τον αριθμό αυτό μπορούμε να είμαστε βέβαιοι ότι τα δεδομένα θέσης τα οποία θα λάβουμε από την διαδικτυακή αυτή βάση δεδομένων θα αναφέρονται στο είδος το οποίο μας ενδιαφέρει.
Στη συνέχεια, θα δούμε πόσα γεωαναφερμένα δεδομένα θέσης έχουμε στην διάθεση μας.
3.3 Μεταφόρτωση των γεωαναφερμένων δεδομένων θέσης
## ===========================================================================##
## Download data
## ===========================================================================##
sp_occ <- occ_search(taxonKey = spp_petromarula$data$key[1], country = "GR", hasCoordinate = T,
limit = 1000)
## ===========================================================================##
## ===========================================================================##
## To prevent any problem with the pathway it is a good practice to remove blank
## spaces from species names
## ===========================================================================##
sp_occ$data$name <- sub(" ", "_", sp_occ$data$name)
## ===========================================================================##
## ===========================================================================##
## Total number of occurrences
## ===========================================================================##
sort(table(sp_occ$data$name), decreasing = T)##
## Petromarula_pinnata (L.) A.DC. Phyteuma_pinnatum L.
## 122 1
Θα πρέπει να τονίσουμε το εξής: Εφ’όσον δεν έχουμε συλλέξει πρωτογενώς τα δεδομένα μας, θα χρειαστεί να ελέγξουμε εάν οι θέσεις εμφάνισης που μεταφορτώσαμε από τις διαδικτυακές βάσεις δεδομένων, συμφωνούν με τα γνωστά όρια εξάπλωσης του είδους που μας ενδιαφέρει. Εάν όντως υπάρχουν κάποια σημεία τα οποία βρίσκονται εκτός της περιοχής εξάπλωσης, θα χρειαστεί να τα αφαιρέσουμε από την ανάλυση.
## ===========================================================================##
## Clean the coordinates
## ===========================================================================##
flags_gbif <- clean_coordinates(x = sp_occ$data, lon = "decimalLongitude", lat = "decimalLatitude",
countries = "countryCode", species = "species", tests = c("capitals", "centroids",
"equal", "gbif", "zeros", "seas"), seas_ref = Crete)
## ===========================================================================##
## ===========================================================================##
## Keep only those that are OK
## ===========================================================================##
flags_gbif <- flags_gbif %>% mutate(Data = ifelse(.summary == "TRUE", "Clean", "Flagged"),
NAME_1 = "Crete", Dataset = "GBIF", longitude = decimalLongitude, latitude = decimalLatitude,
scientific_name = "Petromarula pinnata") %>% dplyr::filter(str_detect(Data, "Cle")) %>%
dplyr::select(scientific_name, longitude, latitude, Dataset)
## ===========================================================================##
## ===========================================================================##
## Create a duplicate for safekeeping
## ===========================================================================##
flags_gbif_dpl <- flags_gbif
## ===========================================================================##Ας φτιάξουμε έναν φάκελο ώστε να αποθηκεύσουμε τα δεδομένα μας και ας τα σώσουμε, προκειμένου να μπορέσουμε να τα φορτώσουμε κάποια άλλη στιγμή:
dir.create("./Species", recursive = TRUE, showWarnings = FALSE)
dir.create("./Species/Petromarula pinnata", recursive = TRUE, showWarnings = FALSE)
saveRDS(flags_gbif, "./Species/ Petromarula pinnata.rds")4 Μεταφόρτωση περιβαλλοντικών μεταβλητών
Τα κλιματικά δεδομένα μας μπορούμε να τα κατεβάσουμε από το WorldClim, ενώ τα υψομετρικά δεδομένα από εδώ και τα δεδομένα σχετικά με την ξηρασία και την εξατμισιοδιαπνοή από εδώ. - Εάν φυσικά δουλεύουμε για μια συγκεκριμένη εργασία (πτυχιακή, μεταπτυχιακή, κτλ) και όχι στα πλαίσια ενός προπτυχιακού μαθήματος επιλογής. Προς το παρόν, θα χρησιμοποιήσουμε κάποιες εντολές για να κατεβάσουμε τα δεδομένα που θέλουμε μέσω της R.
Όπως μπορείτε να διαπιστώσετε, υπάρχουν πάρα πολλά κλιματικά μοντέλα (επί παραδείγματι, CCSM4, HadGEM-2, κτλ) και τέσσερα (4) διαφορετικά κλιματικά σενάρια που έχουν να κάνουν με τη συγκέντρωση των αερίων του θερμοκηπίου στην ατμόσφαιρα εντός του 21ου αιώνα { ενδεικτική βιβλιογραφία ή άλλη ενδεικτική βιβλιογραφία ή ακόμα μια ενδεικτική βιβλιογραφία ή τελευταία ενδεικτική βιβλιογραφία}7.
Συνεπώς, δεν μπορούμε να χρησιμοποιήσουμε μόνο ένα κλιματικό μοντέλο, πόσω μάλλον ένα κλιματικό σενάριο. Στο σημείο αυτό, είναι εύλογο να αρχίσετε να πανικοβάλεστε8. Ευτυχώς όμως, δεν χρειάζεται και ούτε είναι επιστημονικά σωστό, να χρησιμοποιήσουμε όλα τα κλιματικά μοντέλα. Αναλόγως της περιοχής στην οποία εργαζόμαστε, κάποια κλιματικά μοντέλα αναπαριστούν καλύτερα τις περιβαλλοντικές συνθήκες σε σχέση με κάποια άλλα (McSweeney et al., 2015). Καθώς θα διαπιστώσετε (όταν με το καλό διαβάσετε την προαναφερθείσα εργασία), ακόμα και έτσι, είναι πολλά τα κλιματικά μοντέλα για την περιοχή μας. Όμως, μόνο 3 από αυτά έχουν δεδομένα και για τα 4 κλιματικά σενάρια (BCC, HadGEM, CCSM4). Εμείς χάριν συντομίας, θα εργαστούμε όμως μόνο με ένα (1) κλιματικό μοντέλο και ενα (1) κλιματικό σενάριο.
## ============================================================================##
## First, let's download current climate data for the whole planet
## ============================================================================##
current_climate <- getData("worldclim", var = "bio", res = 2.5)
## ============================================================================##
## ============================================================================##
## Then, let's download future climate data
## ============================================================================##
future_climate <- getData("CMIP5", var = "bio", res = 2.5, rcp = 85, model = "CC",
year = 70)
## ============================================================================##Ας αλλάξουμε τα ονόματα των περιβαλλοντικών μεταβλητών.
names(current_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6",
"bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14",
"bio_15", "bio_16", "bio_17", "bio_18", "bio_19")
names(future_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6",
"bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14",
"bio_15", "bio_16", "bio_17", "bio_18", "bio_19")Καθώς τα κλιματικά μας δεδομένα όπως διαπιστώσατε έχουν εξαιρετικά μεγάλο εύρος τιμών, θα χρειαστεί να τα μετατρέψουμε σε oC όσον αφορά την θερμοκρασία (bio_1). Εδώ μπορείτε να διαπιστώσετε το γιατί9. Αυτό γίνεται με την ακόλουθη εντολή:
gain(current_climate$bio_1) = 0.1
gain(future_climate$bio_1) = 0.1Στη συνέχεια, θα κατεβάσουμε τα υψομετρικά δεδομένα για την χώρα που μας ενδιαφέρει:
altitude <- getData("alt", country = "GRC", mask = TRUE)4.1 Κοπή των περιβαλλοντικών δεδομένων στην έκταση της Ελλάδας
Στη συνέχεια, θα ‘κόψουμε’ τις περιβαλλοντικές μας μεταβλητές στο σχήμα της Ελλάδας και ακολούθως θα δημιουργήσουμε ένα νέο αρχείο που θα περιέχει τις μεταβλητές αυτές και την περιοχή μελέτης μας:
## ============================================================================##
## First for Greece
## ============================================================================##
current_climate_crop <- crop(current_climate, Greece, snap = "in") %>% mask(., Greece)
altitude <- crop(altitude, Greece, snap = "in") %>% mask(., Greece) %>%
## We are forcing the extent of altitude to be exactly the same as that of the
## climate data
resample(., current_climate_crop, "bilinear")
future_climate <- crop(future_climate, Greece, snap = "in") %>% mask(., Greece)
## ============================================================================##
## ============================================================================##
## Merge (stack) climate and altitude
## ============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude)
future_climate_crop <- stack(future_climate, altitude)
## ============================================================================##
## ============================================================================##
## Then for Crete
## ============================================================================##
current_climate_Crete <- crop(current_climate_crop, Crete, snap = "in") %>% mask(.,
Crete)
future_climate_Crete <- crop(future_climate_crop, Crete, snap = "in") %>% mask(.,
Crete)
## ============================================================================##4.2 Γραφική απεικόνιση
Μπορούμε να αναπαραστήσουμε γραφικά τα αποτελέσματα ως εξής:
gplot(current_climate_Crete[[1]]) + geom_raster(aes(fill = value)) + scale_fill_gradientn(colours = c("brown",
"red", "yellow", "darkgreen", "green")) + coord_equal() + xlab("Longitude") +
ylab("Latitude") + ggtitle("Temperature")histogram(current_climate_Crete[[1]])levelplot(current_climate_Crete[[1]])Κάντε το ίδιο για την Αττική, τον Άθω και την Θεσσαλία.
4.3 Βασικά στατιστικά στοιχεία
Στη συνέχεια, θα υπολογίσουμε κάποια βασικά στατιστικά στοιχεία για κάθε αρχείο το οποίο δημιουργήσαμε:
cellStats(current_climate_crop, "quantile")## bio_1 bio_2 bio_3 bio_4 bio_5 bio_6 bio_7 bio_8 bio_9 bio_10 bio_11
## 0% 2.4 55 27 4730.00 180 -95 182 -37 94 106 -53
## 25% 12.2 92 33 5988.75 280 -11 262 60 205 207 37
## 50% 14.5 101 34 6516.00 298 12 288 81 229 230 60
## bio_12 bio_13 bio_14 bio_15 bio_16 bio_17 bio_18 bio_19 GRC_msk_alt
## 0% 385 52.00 0 20 138.0 2 2 110 -2.484375
## 25% 542 74.00 8 32 206.0 34 39 184 111.520996
## 50% 653 103.00 15 50 269.5 65 69 250 289.577637
## [ reached getOption("max.print") -- omitted 2 rows ]
cellStats(current_climate_crop, "mean")## bio_1 bio_2 bio_3 bio_4 bio_5 bio_6
## 14.03543 97.27204 34.09507 6455.04935 295.72570 15.26490
## bio_7 bio_8 bio_9 bio_10 bio_11 bio_12
## 280.46080 80.48347 222.14066 224.08342 60.01553 695.86060
## bio_13 bio_14 bio_15 bio_16 bio_17 bio_18
## 110.14041 15.94489 51.10471 296.70353 61.12325 66.01879
## bio_19 GRC_msk_alt
## 269.81864 430.21763
cellStats(current_climate_crop, "sd")## bio_1 bio_2 bio_3 bio_4 bio_5 bio_6
## 2.854903 14.953890 2.226295 671.676911 25.702217 37.634237
## bio_7 bio_8 bio_9 bio_10 bio_11 bio_12
## 35.202266 31.311853 27.093454 26.693768 33.891841 189.062736
## bio_13 bio_14 bio_15 bio_16 bio_17 bio_18
## 41.171401 9.894798 20.722858 109.157622 32.687959 35.129005
## bio_19 GRC_msk_alt
## 102.446221 401.097651
Κάντε το ίδιο για την Αττική, τον Άθω και την Θεσσαλία.
Μπορούμε να δημιουργήσουμε υποσύνολα των περιοχών που μας ενδιαφέρουν, ώστε να ανταποκρίνονται σε συγκεκριμένες κλιματικές (ή άλλες συνθήκες):
subset_GR <- current_climate_crop[[1]] < 15 & current_climate_crop[[1]] > 10.1
plot(subset_GR)Δοκιμάστε το και για άλλες περιοχές και μεταβλητές.
4.4 Σύγκριση
Τώρα, θα πρέπει να ελέγξουμε εάν τα υψομετρικά δεδομένα έχουν την ίδια ανάλυση (resolution), μέγεθος (nrow, ncol) και έκταση (bbox), με τα περιβαλλοντικά μας δεδομένα:
altitude## class : RasterLayer
## dimensions : 165, 246, 40590 (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : 19.375, 29.625, 34.83333, 41.70833 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : GRC_msk_alt
## values : -2.484375, 2346.14 (min, max)
current_climate_crop## class : RasterStack
## dimensions : 165, 246, 40590, 20 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : 19.375, 29.625, 34.83333, 41.70833 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : bio_1, bio_2, bio_3, bio_4, bio_5, bio_6, bio_7, bio_8, bio_9, bio_10, bio_11, bio_12, bio_13, bio_14, bio_15, ...
## min values : 2.400000, 55.000000, 27.000000, 4730.000000, 180.000000, -95.000000, 182.000000, -37.000000, 94.000000, 106.000000, -53.000000, 385.000000, 52.000000, 0.000000, 20.000000, ...
## max values : 19.30, 123.00, 41.00, 7606.00, 348.00, 104.00, 341.00, 155.00, 270.00, 270.00, 133.00, 1336.00, 226.00, 49.00, 102.00, ...
4.5 Resolution change
Στην περίπτωση κατά την οποία τα δύο σετ δεδομένων έχουν διαφορετική ανάλυση, έκταση και μέγεθος, θα χρειαστεί να αλλάξουμε την ανάλυση των υψομετρικών δεδομένων10.
agg_altitude <- aggregate(altitude,
5, ## Why choose 5 and not any other number?
fun = mean) 4.6 Pixel size
Μετά, θα αλλάξουμε το μέγεθος των pixel των υψομετρικών δεδομένων για να συμφωνεί με αυτό των περιβαλλοντικών δεδομένων:
altitude_rsmp <- resample(agg_altitude, current_climate_crop, method = "bilinear")4.7 Final stack
Τέλος, θα τα ενώσουμε σε ένα αρχείο, τόσο για το παρόν, όσο και για το μέλλον
## ============================================================================##
## Merge (stack) climate and altitude
## ============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude_rsmp)
future_climate_crop <- stack(future_climate_crop, altitude_rsmp)
## ============================================================================##4.8 Αποθήκευση σε συμπιεσμένη μορφή
Καλό είναι να αποθηκεύσετε τα αρχεία αυτά υπό την μορφή Rds και να τα αφαιρέσετε από την μνήμη της R11. Αυτό το κάνουμε, προκειμένου να μπορούμε εύκολα να το φορτώνουμε στην R και να μην υπάρχουν τα δεδομένα αυτά συνεχώς στην μνήμη της session που έχουμε ανοικτή, διότι αυξάνουν αρκετά τον χρόνο απόκρισης της (καθυστερεί πολύ η R). Μπορείτε ανά πάσα στιγμή να τα φορτώσετε στην R, με την ακόλουθη εντολή:
## Creating a folder to save our data
dir.create("./Data/", recursive = TRUE, showWarnings = FALSE)
## Save our data in rds format
saveRDS(current_climate_crop, "./Data/Environmental variables and Elevation.rds")
saveRDS(future_climate_crop, "./Data/Environmental variables and Elevation 2070.rds")
## Load our data from rds format
var_GR <- readRDS("./Data/Environmental variables and Elevation.rds")
var_GR_70 <- readRDS("./Data/Environmental variables and Elevation 2070.rds")4.9 Έλεγχος συγγραμικότητας και Variance Inflation Factor (VIF)
Όπως θα γνωρίζετε, ΔΕΝ μπορούμε να χρησιμοποιήσουμε σε μια ανάλυση δύο μεταβλητές, οι οποίες εμφανίζουν συντελεστή συσχέτισης \(> 0.7\)11, καθότι στην προκειμένη περίπτωση οι δύο αυτές μεταβλητές είναι συγγραμικές και ως εκ τούτου εάν τις συμπεριλάβουμε και τις δύο στην ανάλυση μας, τότε δεν θα γνωρίζουμε σε ποιά εκ των δύο θα οφείλεται το αποτέλεσμα της ανάλυσης μας.
Προκειμένου να ξεπεράσουμε αυτόν τον σκόπελο, θα τρέξουμε την ακόλουθη εντολή12:
.
uncorrelated_vars <- vifcor(current_climate_Crete %>% as.data.frame(), th = 0.7)
uncorrelated_vars <- uncorrelated_vars@results$Variables[uncorrelated_vars@results$VIF <=
5]
current_climate_Crete <- stack(subset(current_climate_Crete, uncorrelated_vars))
future_climate_Crete <- stack(subset(future_climate_Crete, uncorrelated_vars))#Δεδομένα θέσης και περιβαλλοντικές μεταβλητές
Πλέον μπορούμε να συνδέσουμε τις γεωαναφερμένες θέσεις εμφάνισης του είδους που μας ενδιαφέρει με δεδομένα ψηφιδοπλέγματος. Στο σημείο αυτό, μπορούμε να εξάγουμε τα περιβαλλοντικά δεδομένα για κάθε μια εκ των θέσεων εμφάνισης του είδους που μας ενδιαφέρει με μια μόνο εντολή, αναδεικνύοντας κατ’αυτόν τον τρόπο το πλεονέκτημα του να έχουμε ενώσει προηγουμένως τα raster αρχεία μας σε ένα RasterStack.
flags_gbif <- flags_gbif %>% as.data.frame()
coordinates(flags_gbif) <- c("longitude", "latitude")
pts.clim <- raster::extract(current_climate_Crete, flags_gbif, method = "bilinear")
sp_occ_clim <- data.frame(cbind(coordinates(flags_gbif), pts.clim, flags_gbif@data))Και τώρα, ήρθε η στιγμή να οπτικοποιήσουμε τα αποτελέσματα μας και δη, την παρουσία του είδους μας στην περιοχή μελέτης:
map.ext <- extent(Crete)
plot(current_climate_crop[[20]], col = terrain.colors(100), ext = map.ext)
plot(flags_gbif, pch = 16, cex = 0.5, add = T)
plot(Crete, add = T)Στο σημείο αυτό, θα ήθελα να σας τονίσω ότι ορισμένες φορές, πέρα από τους ελέγχους που διενεργήσαμε έως τώρα, ίσως χρειαστεί να προχωρήσουμε σε μια αραίωση (thinning) των θέσεων εμφάνισης. Αυτό μπορεί να είναι απαραίτητο όταν έχουμε πάρα πολλά σημεία τα οποία είναι αρκετά κοντά το ένα στο άλλο13. Σε κάθε περίπτωση, το επιστημονικά ορθό είναι να υπάρχει μια (1) θέση εμφάνισης ανά km2. Ένα πακέτο που μπορεί να μας βοηθήσει είναι το spThin14.
## Min. distance: 1000 m Number of reps: 1000
spoccs <- flags_gbif_dpl %>% as.data.frame() %>% dplyr::rename(x = longitude, y = latitude,
Species = scientific_name) %>% addmainfields(species = "Species")
thinned <- thin(loc.data = spoccs, lat.col = "y", long.col = "x", spec.col = "Species",
thin.par = 5, reps = 1000, verbose = T, out.dir = "./Thinned", locs.thinned.list.return = TRUE,
write.files = TRUE, out.base = "Petromarula pinnata thinned")## **********************************************
## Beginning Spatial Thinning.
## Script Started at: Sun Jun 06 18:46:29 2021
## lat.long.thin.count
## 52
## 1000
## [1] "Maximum number of records after thinning: 52"
## [1] "Number of data.frames with max records: 1000"
## [1] "Writing new *.csv files"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin1.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin2.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin3.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin4.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin5.csv"
5 Μοντελοποίηση
Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο εδώ
5.1 Εισαγωγή
Είμαστε πλεόν στην ευχάριστη θέση να έχουμε στην διάθεση μας όλα τα δεδομένα τα οποία είναι απαραίτητα για να προχωρήσουμε στο κυρίως σκέλος της ανάλυσης μας. Κανονικά θα έπρεπε να χρησιμοποιήσουμε πολλά μοντέλα, τα οποία ανήκουν στην ίδια ομάδα (Παρουσίας/Απουσίας ή Ψεύδο-απουσίας - Presence/Absence: P/A), αλλά και σε άλλη ομάδα (Μόνο Παρουσίας - Presence Only: P/O). Κάθε ένα εκ των μοντέλων αυτών έχει τα πλεονεκτήματα και τα μειονεκτήματα του15. Στο σημείο αυτό, πρέπει να τονίσουμε ότι είναι εντελώς λάθος να χρησιμοποιούμε δύο μοντέλα16 τα οποία ανήκουν σε διαφορετικές ομάδες (εν αντιθέσει με ότι προτείνεται στο παράδειγμα του προαναφερθέντος βιβλίου), αλλά να χρησιμοποιεί πολλά μοντέλα τα οποία ανήκουν στην ίδια ομάδα (π.χ., P/A), προκειμένου να μπορεί να βγάλει ασφαλή συμπεράσματα σχετικά με το ερώτημα το οποίο καλείται να απάντησει και να είναι σίγουρος ότι δεν πρόκειται περί στατιστικών κατασκευασμάτων (Araújo & New, 2007).
5.2 Φόρτωση των δεδομένων θέσης
Τώρα θα φορτώσουμε ένα από τα set δεδομένων17 τα οποία έχουν προέλθει από το thinning των θέσεων εμφάνισης18.
thinned_1 <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Fourth_Lab/Fourth_Lab/Thinned/Petromarula pinnata thinned_thin1.xlsx")5.3 Μορφοποιήση των δεδομένων
Το πρώτο πράγμα το οποίο θα χρειαστεί να κάνουμε, είναι να μορφοποιήσουμε με τον κατάλληλο τρόπο τα δεδομένα μας, ώστε να μπορέσουμε να τρέξουμε τις εντολές από τα πακέτα biomod και ecospat. Το δεύτερο πακέτο θα το χρησιμοποιήσουμε στην περίπτωση όπου οι θέσεις εμφάνισης του είδους που μας ενδιαφέρει είναι \(<50\)19.
Ένα άλλο σημείο το οποιο πρέπει να προσέξουμε, έχει να κάνει με τον αριθμό των απουσιών (absences - στην πλειονότητα των περιπτώσεων βέβαια πρόκειται για ψευδοαπουσίες -pseudoabsences, καθώς στην πραγματικότητα δεν έχουμε δεδομένα απουσίας του εκάστοτε είδους). Ο αριθμός αυτός πρέπει να είναι ίσος με τον αριθμό των παρουσιών20.
Τέλος, ένα τρίτο, κρίσιμο σημείο είναι ο αριθμός των set των pseudoabsences που θα περιλαμβάνει το μοντέλο μας, προκειμένου να αποφύγουμε τον σκόπελο της μεροληπτικής δειγματοληψίας (sampling bias)21.
sp_formated_data <- BIOMOD_FormatingData(resp.var = rep(1, nrow( thinned_1 ) ),
## Our rasterStack with the uncorrelated
## variables
expl.var = current_climate_Crete,
## The names of our coords
resp.xy = thinned_1[,c('x', 'y')],
## Species' name
resp.name = "Petromarula pinnata",
## Number of pseudoabsences' sets
PA.nb.rep = 2,
## Number equaling that of the presences
PA.nb.absences = nrow(thinned_1),
## You can change this to 'disk' if you'd
## like
PA.strategy = 'random',
## Min distance from occurrences
PA.dist.min = NULL,
## Max distance from occurrences
PA.dist.max = NULL) ##
## -=-=-=-=-=-=-=-=-=-=-= Petromarula pinnata Data Formating -=-=-=-=-=-=-=-=-=-=-=
##
## Response variable name was converted into Petromarula.pinnata
## ! No data has been set aside for modeling evaluation
## > Pseudo Absences Selection checkings...
## > random pseudo absences selection
## > Pseudo absences are selected in explanatory variables
## ! Some NAs have been automatically removed from your data
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
sp_formated_data##
## -=-=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.formated.data.PA' -=-=-=-=-=-=-=-=-=-=-=-=-=
##
## sp.name = Petromarula.pinnata
##
## 31 presences, 0 true absences and 75 undifined points in dataset
##
##
## 4 explanatory variables
##
## bio_3 bio_5 bio_13 bio_15
## Min. :29.00 Min. :222.0 Min. :111.0 Min. :80.00
## 1st Qu.:30.00 1st Qu.:269.0 1st Qu.:138.2 1st Qu.:83.00
## Median :31.00 Median :277.0 Median :153.0 Median :84.00
## Mean :31.24 Mean :275.0 Mean :160.3 Mean :84.02
## 3rd Qu.:32.00 3rd Qu.:286.8 3rd Qu.:184.8 3rd Qu.:85.00
## Max. :33.00 Max. :301.0 Max. :213.0 Max. :89.00
##
##
## 2 Pseudo Absences dataset available ( PA1 PA2 ) with 39
## absences in each (true abs + pseudo abs)
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Ας αναπαραστήσουμε γραφικά ότι φτιάξαμε στο προηγούμενο βήμα:
plot(sp_formated_data)Η εικόνα αυτή μας δείχνει που βρίσκονται οι ψευδοαπουσίες (pseudoabsences) στο set των ψευδοαπουσιών σε σχέση με τις θέσεις εμφάνισης. Όπως αναμενόταν, οι ψευδοαπουσίες είναι τυχαία κατανεμημένες.
5.4 Παραμετροποιήση των μοντέλων
Πλέον μπορούμε να παραμετροποιήσουμε τα μοντέλα που θέλουμε να τρέξουμε:
modelling_options <- BIOMOD_ModelingOptions()Τώρα, προσδιορίζουμε ποιά μοντέλα θα τρέξουμε, καθώς και πολλές άλλες λεπτομέρειες:
## Always change the name of this object, so as to match the name of your species
Petromarula_pinnata_models <- BIOMOD_Modeling( data = sp_formated_data,
## Here, we indicate which models we are
## going to run (remember that they need
## to be from the same group P/A ή P/O).
## Here, we are going to use some models
## belonging to the P/A group
## ('GLM','GBM','GAM','CTA','ANN','SRE',
## 'FDA','MARS','RF')
models = c('RF', 'CTA'),
models.options = modelling_options,
## How many times we will run the
## analysis in order to evaluate its'
## results [min: 10 - usually: 100]
NbRunEval = 1,
## Training-Testing data split
DataSplit = 80,
## How many times we will run the
## analysis in order to evaluate the
## importance of each IV [min: 10]
VarImport = 1,
## Evaluation criteria
models.eval.meth = c('TSS','ROC','KAPPA'),
## Whether or not we will save the results
SaveObj = TRUE,
## Should the models be re-scaled in 0-1000 range?
rescal.all.models = TRUE,
## Should our models be weighted based on all the data?
do.full.models = TRUE,
## The name of our analysis
modeling.id = "allmodels") Ας φτιάξουμε τους φακέλους όπου θα αποθηκεύσουμε τα μοντέλα μας:
## Create some folders to store our results
dir.create("./BIOMOD2 Modelling Outpout/", recursive = TRUE, showWarnings = FALSE)
dir.create("./BIOMOD2 Modelling Outpout/Species/", recursive = TRUE, showWarnings = FALSE)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata", recursive = TRUE,
showWarnings = FALSE)
dir.create("./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current", recursive = TRUE,
showWarnings = FALSE)
## Save the results in rds format
saveRDS(Petromarula_pinnata_models, "./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current/ Petromarula_pinnata models.rds")5.5 Αξιολόγηση των αρχικών μοντέλων
Ας δημιουργήσουμε ένα αντικείμενο το οποίο θα περιέχει την αξιολόγηση των μοντέλων:
Petromarula_pinnata_scores <- get_evaluations(Petromarula_pinnata_models)5.6 Σημαντικότητα μεταβλητών
Ας δούμε τη σημαντικότητα των μεταβλητών για κάθε μοντέλο (RF, CTA), κάθε σετ δεδομένων (PA1-2) και κάθε επανάληψη (Run, Full):
(Petromarula_pinnata_var_import <- get_variables_importance(Petromarula_pinnata_models))Ας δούμε τη σημαντικότητα των μεταβλητών για κάθε μοντέλο (RF, CTA) συνολικά (δηλαδή τον μέσο όρο):
options(digits = 2)
apply(Petromarula_pinnata_var_import, c(1, 2), mean)## RF CTA
## bio_3 0.17 0.355
## bio_5 0.33 0.431
## bio_13 0.37 0.632
## bio_15 0.17 0.086
5.7 Ensemble modelling
Τώρα θα τρέξουμε την ανάλυση που θα ενώνει τα αποτελέσματα των δύο μοντέλων (ονομάζεται ensemble modelling αυτή η τεχνική ανάλυσης):
Petromarula_pinnata_ensemble_models <- BIOMOD_EnsembleModeling( modeling.output = Petromarula_pinnata_models,
em.by = 'all',
chosen.models = 'all',
# You can change that
eval.metric = 'TSS',
# Keep only the models with
# TSS >= 0.8
eval.metric.quality.threshold = 0.8,
models.eval.meth = c('KAPPA','TSS','ROC'),
prob.mean = T,
committee.averaging = TRUE,
prob.mean.weight = T,
prob.mean.weight.decay = 'proportional',
VarImport = 0 )dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble",
recursive = TRUE, showWarnings = FALSE)
saveRDS(Petromarula_pinnata_ensemble_models, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble/Petromarula_pinnata ensemble.rds")5.8 Αξιολόγηση των συνολικών μοντέλων (ensemble modelling)
Ας δημιουργήσουμε ένα αντικείμενο το οποίο θα περιέχει την αξιολόγηση των μοντέλων από το ensemble modelling. Όσο πιο κοντά στη μονάδα είναι οι τιμές για τους δείκτες TSS, ROC και KAPPA, τόσο καλύτερα22:
options(digits = 4)
(Petromarula_pinnata_ensemble_models_scores <- get_evaluations(Petromarula_pinnata_ensemble_models))## $Petromarula.pinnata_EMmeanByTSS_mergedAlgo_mergedRun_mergedData
## Testing.data Cutoff Sensitivity Specificity
## KAPPA 0.911 616.0 96.77 96
## TSS 0.928 616.0 96.77 96
## ROC 0.994 618.5 96.77 96
##
## $Petromarula.pinnata_EMcaByTSS_mergedAlgo_mergedRun_mergedData
## Testing.data Cutoff Sensitivity Specificity
## KAPPA 0.909 747 93.55 97.33
## TSS 0.909 747 93.55 97.33
## ROC 0.969 750 93.55 97.33
##
## $Petromarula.pinnata_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData
## Testing.data Cutoff Sensitivity Specificity
## KAPPA 0.911 611.0 96.77 96
## TSS 0.928 611.0 96.77 96
## ROC 0.994 611.5 96.77 96
## ===========================================================## Vizualisation
## ===========================================================##
g2 <- models_scores_graph(Petromarula_pinnata_ensemble_models, by = "models", metrics = c("TSS",
"KAPPA"))5.9 Προβολή στο μέλλον
Ας τρέξουμε την ανάλυση για το πώς θα αλλάξει η κατάσταση του είδους που μας ενδιαφέρει στο μέλλον. Αρχικά θα τρέξουμε την ανάλυση για κάθε μοντέλο ξεχωριστά (BIOMOD_Projection) και συνολικά (BIOMOD_EnsembleForecasting) για τις παρούσες συνθήκες και στη συνέχεια για τις μελλοντικές συνθήκες (με τις ίδιες εντολές, αλλά θα αλλάξουμε το raster αρχείο στην εντολή new.env):
Petromarula_pinnata_models_proj_current <- BIOMOD_Projection(modeling.output = Petromarula_pinnata_models,
new.env = current_climate_Crete, proj.name = "current", binary.meth = "TSS",
selected.models = "all", build.clamping.mask = TRUE, output.format = ".img",
do.stack = T)Petromarula_pinnata_ensemble_models_proj_current <- BIOMOD_EnsembleForecasting(EM.output = Petromarula_pinnata_ensemble_models,
projection.output = Petromarula_pinnata_models_proj_current, binary.meth = "TSS",
output.format = ".img", do.stack = T)dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/",
recursive = TRUE, showWarnings = FALSE)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual",
recursive = TRUE, showWarnings = FALSE)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble",
recursive = TRUE, showWarnings = FALSE)
saveRDS(Petromarula_pinnata_models_proj_current, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual/Petromarula_pinnata Current Individual Projections.rds")
saveRDS(Petromarula_pinnata_ensemble_models_proj_current, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble/Petromarula_pinnata Current Ensemble Projections.rds")Petromarula_pinnata_proj_2070_cc85 <- BIOMOD_Projection(modeling.output = Petromarula_pinnata_models,
new.env = future_climate_Crete, proj.name = "2070_CC85", binary.meth = "TSS",
output.format = ".img", do.stack = T)dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070", recursive = TRUE,
showWarnings = FALSE)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/",
recursive = TRUE, showWarnings = FALSE)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual",
recursive = TRUE, showWarnings = FALSE)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble",
recursive = TRUE, showWarnings = FALSE)
saveRDS(Petromarula_pinnata_proj_2070_cc85, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual/Petromarula_pinnata 2070 Individual Projections.rds")Petromarula_pinnata_ensemble_models_proj_2070_cc85 <- BIOMOD_EnsembleForecasting(EM.output = Petromarula_pinnata_ensemble_models,
projection.output = Petromarula_pinnata_proj_2070_cc85, binary.meth = "TSS",
output.format = ".img", do.stack = T)saveRDS(Petromarula_pinnata_ensemble_models_proj_2070_cc85, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble/Petromarula_pinnata 2070 Ensemble Projections.rds")5.10 Οπτικοποιήση των συνολικών μοντέλων
Τέλος, ας οπτικοποιήσουμε το αποτέλεσμα της ανάλυσης μας. Πρώτα θα θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει την αξιολόγηση της κατάστασης της Petromarula pinnata μέσω της τεχνικής ensemble modelling. Στη συνέχεια, θα κρατήσουμε μόνο τα raster layers για το mean και το weighted mean. Τέλος, θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει μόνο τις συντεταγμένες των θέσεων εμφάνισης του είδους που μας ενδιαφέρει:
stk_Petromarula_pinnata_ef_2070_cc85 <- get_predictions(Petromarula_pinnata_ensemble_models_proj_2070_cc85)
stk_Petromarula_pinnata_ef_2070_cc85 <- subset(stk_Petromarula_pinnata_ef_2070_cc85,
grep("EMmean|EMwmean", names(stk_Petromarula_pinnata_ef_2070_cc85)))
names(stk_Petromarula_pinnata_ef_2070_cc85) <- sapply(strsplit(names(stk_Petromarula_pinnata_ef_2070_cc85),
"_"), getElement, 2)
Petromarula_pinnata_coords <- thinned_1 %>% dplyr::select(x, y)Στη συνέχεια θα οπτικοποιήσουμε το αποτέλεσμα της ανάλυσης μας έως αυτό το σημείο:
library(viridis)
levelplot(stk_Petromarula_pinnata_ef_2070_cc85, main = "Petromarula pinnata ensemble projections\nin 2070 with cc85",
col.regions = viridis(256)) + layer(sp.points(SpatialPoints(Petromarula_pinnata_coords),
pch = 20, col = 1))5.11 Μεταβολή περιοχής εξάπλωσης
Έχουμε φτάσει πλέον στο τελικό σημείο της ανάλυσης μας, όπου θα φορτώσουμε τα αποτελέσματα της ανάλυσης για το mean και το weighted mean της Petromarula pinnata και στη συνέχεια θα υπολογίσουμε και θα οπτικοποιήσουμε τη μεταβολή της περιοχής εξάπλωσης της στο μέλλον:
##==================================
## Load the the data
##==================================
## Current
Petromarula_pinnata_bin_proj_current <- stack('Petromarula.pinnata/proj_current/proj_current_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_current) <- c('ca', 'm', 'wm')
## 2070
Petromarula_pinnata_bin_proj_2070_CC85 <- stack('Petromarula.pinnata/proj_2070_CC85/proj_2070_CC85_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_2070_CC85) <- c('ca', 'm', 'wm')
##=================================
## Species Range Change
##=================================
SRC_current_2070_CC85 <- BIOMOD_RangeSize( Petromarula_pinnata_bin_proj_current,
Petromarula_pinnata_bin_proj_2070_CC85 )
SRC_current_2070_CC85$Compt.By.Models## Loss Stable0 Stable1 Gain PercLoss PercGain SpeciesRangeChange
## ca 103 303 27 33 79.23 25.39 -53.85
## m 106 330 12 18 89.83 15.25 -74.58
## wm 100 297 31 38 76.34 29.01 -47.33
## CurrentRangeSize FutureRangeSize.NoDisp FutureRangeSize.FullDisp
## ca 130 27 60
## m 118 12 30
## wm 131 31 69
Petromarula_pinnata_src_map <- SRC_current_2070_CC85$Diff.By.Pixel
names(Petromarula_pinnata_src_map) <- c("ca cur-2070", "m cur-2070", "wm cur-2070")
##=================================
## Vizualisation
##=================================
my.at <- seq(-2.5,1.5,1)
myColorkey <- list(
## color change
at = my.at,
labels = list(
## labels
labels = c("lost", "pres", "abs","gain"),
## legend
at = my.at[-1]-0.5
))
rasterVis::levelplot( SRC_current_2070_CC85$Diff.By.Pixel,
main = "Petromarula pinnata range change",
colorkey = myColorkey,
layout = c(1,3),
)6 Εργασία για το σπίτι
Έχετε διορία δεκατεσσάρων (14) ημερών να στείλετε την εργασία σας σε μορφή PDF23 σε αυτό το e-mail.
1. Κατεβάστε γεωγραφικά δεδομένα για την Ελλάδα.
2. Δημιουργήστε ξεχωριστά αρχεία για την Κρήτη, την Αττική, τον Άθω, την Μακεδονία και τη Θεσσαλία.
3. Κατεβάστε τα απαραίτητα κλιματικά δεδομένα.
4. Κατεβάστε υψομετρικά δεδομένα για την Ελλάδα.
5. Υπολογίστε τις βασικές στατιστικές παραμέτρους για όλα τα αρχεία τα οποία δημιουργήσατε.
6. Υπολογίστε τον συντελεστή συσχέτισης για όλα τα περιβαλλοντικά δεδομένα για όλα τα αρχεία που έχετε δημιουργήσει. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.7\), θα τις χρησιμοποιήσετε όλες;
7. Κατεβάστε δεδομένα θέσης για τα είδη Asperula pubescens, Centaurea idaea και Scutellaria sieberi.
8. Παραθέστε πληροφορίες σχετικά με τα είδη αυτά (Καθεστώς κινδύνου/προστασίας, καθεστώς ενδημισμού, ταξινομικές πληροφορίες). Προχωρήστε στις απαραίτητες ενέργειες, έτσι ώστε να διαπιστώσετε εάν τα δεδομένα αυτά εντοπίζονται στη σωστή γεωγραφικά περιοχή.
9. Πόσες θέσεις εμφάνισης αριθμούσαν τα είδη αυτά πριν και μετά την διαδικασία αραίωσης;
10. Φορτώστε τα τελικά δεδομένα θέσης για κάθε είδος για το οποίο έχετε κατεβάσει δεδομένα θέσης.
11. Αλλάξτε τα rasterStack αρχεία με αυτά για την Κρήτη σε όλα τα απαραίτητα σημεία του κώδικα. Αλλάξτε επίσης το όνομα του είδους σας. Αλλάξτε τον αριθμό των σετ σε 10 και τον αριθμό των ψευδοαπουσιών24. Αλλάξτε τον αριθμο των φορών που θα τρέξει η ανάλυση σε 10.
12. Ποιά είναι η σημαντικότερη μεταβλητή σε κάθε αλγόριθμο; Συμπίπτουν μεταξύ των αλγορίθμων;
13. Ποιό είναι το καλύτερο μοντέλο σύμφωνα με τους δείκτες αξιολόγησης TSS και KAPPA;.
14. Προβάλετε τα μοντέλα σας στο μέλλον. Παραθέστε τους σχετικούς χάρτες.
15. Πώς μεταβάλλεται η περιοχή εξάπλωσης των ειδών αυτών στο μέλλον;25 Τι συμβαίνει στην περίπτωση της μηδενικής και της πλήρους διασποράς, αντίστοιχα;
16. Παραθέστε το σχετικό γράφημα.
7 Appendix: R code
##===========================================================================##
## Do NOT run this code
##===========================================================================##
library(knitr)
library(rmdformats)
## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
cache=TRUE,
prompt=FALSE,
tidy=TRUE,
# comment=NA,
message=FALSE,
warning=FALSE,
eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
pre {
max-height: 400px;
overflow-y: auto;
}
pre[class] {
max-height: 400px;
}
##===========================================================================##
## Install the main packages
##===========================================================================##
install.packages(c('tidyverse', 'rgbif', 'raster', 'data.table',
'rgeos', 'rgdal', 'gridExtra',
'dismo', 'biomod2', 'sf', 'usdm',
'speciesgeocodeR', 'pacman', 'spThin',
'CoordinateCleaner', 'ggspatial', 'scales',
'scrubr', 'mapr', 'rasterVis', 'ade4',
'viridis', 'biogeo'), dependencies = T)
##===========================================================================##
##===========================================================================##
## Load them
##===========================================================================##
pacman::p_load(tidyverse, rgbif, raster, usdm,
rgeos, rgdal, usdm, biomod2, gridExtra,
countrycode, devtools, sf, spThin, biogeo,
speciesgeocodeR, pacman, dismo,
CoordinateCleaner, ggspatial, scales,
scrubr, mapr, rasterVis, data.table, ade4, viridis)
##===========================================================================##
##===========================================================================##
## Download the data for Greece
##===========================================================================##
## There are three levels available, you are free to explore them
Greece <- getData('GADM', country = 'GRC', level = 3)
##===========================================================================##
##===========================================================================##
## Plot Greece
##===========================================================================##
plot(Greece)
##===========================================================================##
## First we need to see what is inside the object we created
##===========================================================================##
Greece
names(Greece)
Crete <- subset(Greece, NAME_1 == "Crete")
##===========================================================================##
## Plot Crete
##===========================================================================##
plot(Crete)
spp_petromarula <- name_suggest(q = 'Petromarula pinnata', rank='species',limit = 10000)
spp_petromarula
##===========================================================================##
## Download data
##===========================================================================##
sp_occ <- occ_search(taxonKey = spp_petromarula$data$key[1],
country ='GR',
hasCoordinate = T,
limit = 1000)
##===========================================================================##
##===========================================================================##
## To prevent any problem with the pathway it is a good practice to remove blank
## spaces from species names
##===========================================================================##
sp_occ$data$name <- sub(" ", "_", sp_occ$data$name)
##===========================================================================##
##===========================================================================##
## Total number of occurrences
##===========================================================================##
sort(table(sp_occ$data$name), decreasing = T)
##===========================================================================##
## Clean the coordinates
##===========================================================================##
flags_gbif <- clean_coordinates(x = sp_occ$data,
lon = "decimalLongitude",
lat = "decimalLatitude",
countries = "countryCode",
species = "species",
tests = c("capitals",
"centroids",
"equal",
"gbif",
"zeros",
"seas"),
seas_ref = Crete)
##===========================================================================##
##===========================================================================##
## Keep only those that are OK
##===========================================================================##
flags_gbif <- flags_gbif %>%
mutate(Data = ifelse(.summary == 'TRUE',
'Clean',
'Flagged'),
NAME_1 = 'Crete',
Dataset = 'GBIF',
longitude = decimalLongitude,
latitude = decimalLatitude,
scientific_name = 'Petromarula pinnata') %>%
dplyr::filter(str_detect(Data, 'Cle')) %>%
dplyr::select(scientific_name,
longitude,
latitude,
Dataset)
##===========================================================================##
##===========================================================================##
## Create a duplicate for safekeeping
##===========================================================================##
flags_gbif_dpl <- flags_gbif
##===========================================================================##
dir.create("./Species", recursive = TRUE, showWarnings = FALSE)
dir.create("./Species/Petromarula pinnata", recursive = TRUE, showWarnings = FALSE)
saveRDS(flags_gbif, "./Species/ Petromarula pinnata.rds")
##============================================================================##
## First, let's download current climate data for the whole planet
##============================================================================##
current_climate <- getData('worldclim', var = 'bio', res = 2.5)
##============================================================================##
##============================================================================##
## Then, let's download future climate data
##============================================================================##
future_climate <- getData('CMIP5', var = 'bio', res = 2.5,
rcp = 85, model = 'CC', year = 70)
##============================================================================##
names(current_climate) <- c('bio_1', 'bio_2', 'bio_3', 'bio_4', 'bio_5',
'bio_6', 'bio_7', 'bio_8', 'bio_9', 'bio_10',
'bio_11', 'bio_12', 'bio_13', 'bio_14', 'bio_15',
'bio_16', 'bio_17', 'bio_18', 'bio_19')
names(future_climate) <- c('bio_1', 'bio_2', 'bio_3', 'bio_4', 'bio_5',
'bio_6', 'bio_7', 'bio_8', 'bio_9', 'bio_10',
'bio_11', 'bio_12', 'bio_13', 'bio_14', 'bio_15',
'bio_16', 'bio_17', 'bio_18', 'bio_19')
gain(current_climate$bio_1) = 0.1
gain(future_climate$bio_1) = 0.1
altitude <- getData('alt', country = 'GRC', mask = TRUE)
##============================================================================##
## First for Greece
##============================================================================##
current_climate_crop <- crop(current_climate,
Greece,
snap = 'in') %>%
mask(., Greece)
altitude <- crop(altitude,
Greece,
snap = 'in') %>%
mask(., Greece) %>%
## We are forcing the extent of altitude to be exactly the same as that of the
## climate data
resample(., current_climate_crop, 'bilinear')
future_climate <- crop(future_climate,
Greece,
snap = 'in') %>%
mask(., Greece)
##============================================================================##
##============================================================================##
## Merge (stack) climate and altitude
##============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude)
future_climate_crop <- stack(future_climate, altitude)
##============================================================================##
##============================================================================##
## Then for Crete
##============================================================================##
current_climate_Crete <- crop(current_climate_crop,
Crete,
snap = 'in') %>%
mask(., Crete)
future_climate_Crete <- crop(future_climate_crop,
Crete,
snap = 'in') %>%
mask(., Crete)
##============================================================================##
gplot(current_climate_Crete[[1]]) +
geom_raster(aes(fill = value)) +
scale_fill_gradientn(colours = c("brown",
"red",
"yellow",
"darkgreen",
"green")) +
coord_equal() +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Temperature")
histogram(current_climate_Crete[[1]])
levelplot(current_climate_Crete[[1]])
cellStats(current_climate_crop, "quantile")
cellStats(current_climate_crop, "mean")
cellStats(current_climate_crop, "sd")
subset_GR <- current_climate_crop[[1]] < 15 & current_climate_crop[[1]] > 10.1
plot(subset_GR)
altitude
current_climate_crop
agg_altitude <- aggregate(altitude,
5, ## Why choose 5 and not any other number?
fun = mean)
altitude_rsmp <- resample(agg_altitude, current_climate_crop, method = "bilinear")
## ============================================================================##
## Merge (stack) climate and altitude
## ============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude_rsmp)
future_climate_crop <- stack(future_climate_crop, altitude_rsmp)
## ============================================================================##
## Creating a folder to save our data
dir.create("./Data/", recursive = TRUE, showWarnings = FALSE)
## Save our data in rds format
saveRDS(current_climate_crop, "./Data/Environmental variables and Elevation.rds")
saveRDS(future_climate_crop, "./Data/Environmental variables and Elevation 2070.rds")
## Load our data from rds format
var_GR <- readRDS("./Data/Environmental variables and Elevation.rds")
var_GR_70 <- readRDS("./Data/Environmental variables and Elevation 2070.rds")
uncorrelated_vars <- vifcor(current_climate_Crete %>% as.data.frame(), th = 0.7)
uncorrelated_vars <- uncorrelated_vars @results$Variables[uncorrelated_vars@results$VIF <= 5]
current_climate_Crete <- stack(subset(current_climate_Crete,
uncorrelated_vars))
future_climate_Crete <- stack(subset(future_climate_Crete,
uncorrelated_vars))
flags_gbif <- flags_gbif %>% as.data.frame()
coordinates(flags_gbif) <- c("longitude", "latitude")
pts.clim <- raster::extract(current_climate_Crete, flags_gbif, method = "bilinear")
sp_occ_clim <- data.frame(cbind(coordinates(flags_gbif), pts.clim, flags_gbif @data))
map.ext <- extent(Crete)
plot(current_climate_crop[[20]], col = terrain.colors(100), ext = map.ext)
plot(flags_gbif, pch = 16, cex = 0.5, add = T)
plot(Crete, add = T)
## Min. distance: 1000 m
## Number of reps: 1000
spoccs <- flags_gbif_dpl %>%
as.data.frame() %>%
dplyr::rename(x = longitude,
y = latitude,
Species = scientific_name) %>%
addmainfields(species = "Species")
thinned <- thin(loc.data = spoccs,
lat.col = "y",
long.col = "x",
spec.col = "Species",
thin.par = 5,
reps = 1000,
verbose = T,
out.dir = "./Thinned",
locs.thinned.list.return = TRUE,
write.files = TRUE, out.base = "Petromarula pinnata thinned")
thinned_1 <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Fourth_Lab/Fourth_Lab/Thinned/Petromarula pinnata thinned_thin1.xlsx")
sp_formated_data <- BIOMOD_FormatingData(resp.var = rep(1, nrow( thinned_1 ) ),
## Our rasterStack with the uncorrelated
## variables
expl.var = current_climate_Crete,
## The names of our coords
resp.xy = thinned_1[,c('x', 'y')],
## Species' name
resp.name = "Petromarula pinnata",
## Number of pseudoabsences' sets
PA.nb.rep = 2,
## Number equaling that of the presences
PA.nb.absences = nrow(thinned_1),
## You can change this to 'disk' if you'd
## like
PA.strategy = 'random',
## Min distance from occurrences
PA.dist.min = NULL,
## Max distance from occurrences
PA.dist.max = NULL)
sp_formated_data
plot(sp_formated_data)
modelling_options <- BIOMOD_ModelingOptions()
## Always change the name of this object, so as to match the name of your species
Petromarula_pinnata_models <- BIOMOD_Modeling( data = sp_formated_data,
## Here, we indicate which models we are
## going to run (remember that they need
## to be from the same group P/A ή P/O).
## Here, we are going to use some models
## belonging to the P/A group
## ('GLM','GBM','GAM','CTA','ANN','SRE',
## 'FDA','MARS','RF')
models = c('RF', 'CTA'),
models.options = modelling_options,
## How many times we will run the
## analysis in order to evaluate its'
## results [min: 10 - usually: 100]
NbRunEval = 1,
## Training-Testing data split
DataSplit = 80,
## How many times we will run the
## analysis in order to evaluate the
## importance of each IV [min: 10]
VarImport = 1,
## Evaluation criteria
models.eval.meth = c('TSS','ROC','KAPPA'),
## Whether or not we will save the results
SaveObj = TRUE,
## Should the models be re-scaled in 0-1000 range?
rescal.all.models = TRUE,
## Should our models be weighted based on all the data?
do.full.models = TRUE,
## The name of our analysis
modeling.id = "allmodels")
## Create some folders to store our results
dir.create('./BIOMOD2 Modelling Outpout/', recursive = TRUE, showWarnings = FALSE)
dir.create('./BIOMOD2 Modelling Outpout/Species/', recursive = TRUE, showWarnings = FALSE)
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata', recursive = TRUE, showWarnings = FALSE)
dir.create('./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current', recursive = TRUE, showWarnings = FALSE)
## Save the results in rds format
saveRDS(Petromarula_pinnata_models, './BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current/ Petromarula_pinnata models.rds')
Petromarula_pinnata_scores <- get_evaluations(Petromarula_pinnata_models)
(Petromarula_pinnata_var_import <- get_variables_importance(Petromarula_pinnata_models))
options(digits = 2)
apply(Petromarula_pinnata_var_import, c(1,2), mean)
Petromarula_pinnata_ensemble_models <- BIOMOD_EnsembleModeling( modeling.output = Petromarula_pinnata_models,
em.by = 'all',
chosen.models = 'all',
# You can change that
eval.metric = 'TSS',
# Keep only the models with
# TSS >= 0.8
eval.metric.quality.threshold = 0.8,
models.eval.meth = c('KAPPA','TSS','ROC'),
prob.mean = T,
committee.averaging = TRUE,
prob.mean.weight = T,
prob.mean.weight.decay = 'proportional',
VarImport = 0 )
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble',
recursive=TRUE, showWarnings=FALSE)
saveRDS(Petromarula_pinnata_ensemble_models,
'./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble/Petromarula_pinnata ensemble.rds')
options(digits = 4)
(Petromarula_pinnata_ensemble_models_scores <- get_evaluations(Petromarula_pinnata_ensemble_models))
##===========================================================##
## Vizualisation
##===========================================================##
g2 <- models_scores_graph( Petromarula_pinnata_ensemble_models,
by = 'models',
metrics = c('TSS','KAPPA'))
Petromarula_pinnata_models_proj_current <- BIOMOD_Projection( modeling.output = Petromarula_pinnata_models,
new.env = current_climate_Crete,
proj.name = "current",
binary.meth = "TSS",
selected.models = 'all',
build.clamping.mask = TRUE,
output.format = ".img",
do.stack = T )
Petromarula_pinnata_ensemble_models_proj_current <-
BIOMOD_EnsembleForecasting( EM.output = Petromarula_pinnata_ensemble_models,
projection.output = Petromarula_pinnata_models_proj_current,
binary.meth = "TSS",
output.format = ".img",
do.stack = T)
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/',
recursive=TRUE, showWarnings=FALSE)
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual',
recursive=TRUE, showWarnings=FALSE)
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble',
recursive=TRUE, showWarnings=FALSE)
saveRDS(Petromarula_pinnata_models_proj_current,
'./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual/Petromarula_pinnata Current Individual Projections.rds')
saveRDS(Petromarula_pinnata_ensemble_models_proj_current,
'./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble/Petromarula_pinnata Current Ensemble Projections.rds')
Petromarula_pinnata_proj_2070_cc85 <- BIOMOD_Projection( modeling.output = Petromarula_pinnata_models,
new.env = future_climate_Crete,
proj.name = "2070_CC85",
binary.meth = "TSS",
output.format = ".img",
do.stack = T )
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070',
recursive=TRUE, showWarnings=FALSE)
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/',
recursive=TRUE, showWarnings=FALSE)
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual',
recursive=TRUE, showWarnings=FALSE)
dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble',
recursive=TRUE, showWarnings=FALSE)
saveRDS(Petromarula_pinnata_proj_2070_cc85,
'./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual/Petromarula_pinnata 2070 Individual Projections.rds')
Petromarula_pinnata_ensemble_models_proj_2070_cc85 <-
BIOMOD_EnsembleForecasting( EM.output = Petromarula_pinnata_ensemble_models,
projection.output = Petromarula_pinnata_proj_2070_cc85,
binary.meth = "TSS",
output.format = ".img",
do.stack = T )
saveRDS(Petromarula_pinnata_ensemble_models_proj_2070_cc85,
'./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble/Petromarula_pinnata 2070 Ensemble Projections.rds')
Petromarula_pinnata_ensemble_models_proj_2070_cc85 <- readRDS('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble/Petromarula_pinnata 2070 Ensemble Projections.rds')
stk_Petromarula_pinnata_ef_2070_cc85 <- get_predictions(Petromarula_pinnata_ensemble_models_proj_2070_cc85)
stk_Petromarula_pinnata_ef_2070_cc85 <- subset(stk_Petromarula_pinnata_ef_2070_cc85,
grep('EMmean|EMwmean',
names(stk_Petromarula_pinnata_ef_2070_cc85)))
names(stk_Petromarula_pinnata_ef_2070_cc85) <- sapply(strsplit(names(stk_Petromarula_pinnata_ef_2070_cc85),
'_'),
getElement, 2)
Petromarula_pinnata_coords <- thinned_1 %>% dplyr::select(x, y)
library(viridis)
levelplot(stk_Petromarula_pinnata_ef_2070_cc85,
main = 'Petromarula pinnata ensemble projections\nin 2070 with cc85',
col.regions = viridis(256)) +
layer(sp.points(SpatialPoints(Petromarula_pinnata_coords), pch=20, col=1))
##==================================
## Load the the data
##==================================
## Current
Petromarula_pinnata_bin_proj_current <- stack('Petromarula.pinnata/proj_current/proj_current_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_current) <- c('ca', 'm', 'wm')
## 2070
Petromarula_pinnata_bin_proj_2070_CC85 <- stack('Petromarula.pinnata/proj_2070_CC85/proj_2070_CC85_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_2070_CC85) <- c('ca', 'm', 'wm')
##=================================
## Species Range Change
##=================================
SRC_current_2070_CC85 <- BIOMOD_RangeSize( Petromarula_pinnata_bin_proj_current,
Petromarula_pinnata_bin_proj_2070_CC85 )
SRC_current_2070_CC85$Compt.By.Models
Petromarula_pinnata_src_map <- SRC_current_2070_CC85$Diff.By.Pixel
names(Petromarula_pinnata_src_map) <- c("ca cur-2070", "m cur-2070", "wm cur-2070")
##=================================
## Vizualisation
##=================================
my.at <- seq(-2.5,1.5,1)
myColorkey <- list(
## color change
at = my.at,
labels = list(
## labels
labels = c("lost", "pres", "abs","gain"),
## legend
at = my.at[-1]-0.5
))
rasterVis::levelplot( SRC_current_2070_CC85$Diff.By.Pixel,
main = "Petromarula pinnata range change",
colorkey = myColorkey,
layout = c(1,3),
)
##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##
klippy::klippy()8 Βιβλιογραφία
Υπάρχουν δύο τρόποι για να το κάνετε αυτό γρήγορα:
1.library(easypackages) packages('package_name_1', 'package_name_2')
2.install.packages(c('package_name_1', 'package_name_2'), dependecies = TRUE)↩︎Ή με αυτήν την εντολή εάν θέλουμε μόνο την Ελλάδα:
ISO_countries <- getData("ISO3") %>% as.data.frame %>% filter(NAME=="Greece")↩︎Ένα τυπικό παράδειγμα είναι ο Άτλαντας της Χλωρίδας του Αιγαίου ή ο Άτλαντας της Χλωρίδας της Ευρώπης, όμως στην περίπτωση αυτή υπεισέρχονται άλλα προβλήματα τα οποία σχετίζονται με την ανάλυση των εν λόγω χαρτών (π.χ., ο Άτλαντας της Χλωρίδας της Ευρώπης είναι σε ανάλυση \(50 \times 50 km^2\))↩︎
Δεν είναι τα μοναδικά όμως: υπάρχουν πολλά άλλα, όπως το dismo και το spocc↩︎
Κάθε ένα από αυτά τα πακέτα, έχει τουλάχιστον μια σύντομη περιγραφή (αγγλιστί: vignette), την οποία είναι ΑΠΑΡΑΙΤΗΤΟ να διαβάσετε - ειδικά αυτές του
rgbif↩︎Με την εντολή
?rgbif::name_suggest()μπορείτε να δείτε τα arguments που δέχεται η εντολή αυτή και τι μπορείτε να αλλάξετε - σας συνιστώ εντόνως να το κάνετε, αλλά και για κάθε εντολή που χρησιμοποιείτε↩︎Καλό είναι να τις διαβάσετε, γιατί θα σας βοηθήσει και στην εργασία σας↩︎
Και εγώ στη θέση σας, το ίδιο θα έκανα↩︎
Γιατί;↩︎
Γιατί όμως να μην αλλάξουμε την ανάλυση των περιβαλλοντικών δεδομένων;↩︎
Αυτό είναι το λεγόμενο hard boundary - άλλοι ερευνητές εφαρμόζουν το όριο \(> 0.8\) ή ακόμα και \(> 0.9\), αλλά είναι επιστημονικά ασφαλέστερο να εφαρμόζει κανείς το αυστηρότερο όριο↩︎
Εδώ εφαρμόζουμε το αυστηρότερο όριο - μπορείτε να το αλλάξετε, εάν επιθυμείτε↩︎
Σε απόσταση \(<1000m\) όταν εργαζόμαστε με χάρτες υποβάθρου κλίμακας 1 x 1 km2 (ή ανάλυσης 30 arc-seconds)↩︎
Προτρέξτε στους Guisan et al. (2017) για περισσότερες λεπτομέρειες - εάν χρειάζεστε περισσότερη βιβλιογραφία, απευθυνθείτε σε εμένα (εάν φυσικά το επιθυμείτε)↩︎
Όπως θα διαπιστώσετε παρακάτω, θα εφαρμόσουμε μια τεχνική η οποία ονομάζεται ensemble modelling. Σύμφωνα με αυτή, είναι προτιμότερο και επιστημονικά ορθότερο να μην χρησιμοποιεί κανείς ένα μόνο μοντέλο (π.χ., GLM ή τo MaxEnt)↩︎
Προσέξτε ότι πρόκειται για αρχείο xlsx και όχι csv - αυτό σημαίνει ότι άνοιξα στο excel το αρχείο csv, χρησιμοποίησα την εντολή
data to columnsκαι το έσωσα ως αρχείο xlsx↩︎Δεν έχει σημασία ποιό θα διαλέξουμε - έχουν όλα την ίδια πιθανότητα↩︎
Θα μπορούσε κάποιος να αντιτείνει ότι μπορούμε να χρησιμοποιήσουμε το πακέτο
ecospatόταν \[\frac{\sum{Ανεξάρτητες,\, μη\, συγγραμικές\, μεταβλητές}}{\sum{Θέσεις\, εμφάνισης}} \geq 0.1\]↩︎Στο βιβλίο των Guisan et al. (2017) αναφέρεται λανθασμένα ότι ο αριθμός των pseudoabsences πρέπει να είναι τουλάχιστον 10000 - ο αριθμός αυτός έχει να κάνει με το MaxEnt, το οποίο ανήκει στην ομάδα των P/O μοντέλων↩︎
Καλό είναι να είναι τουλάχιστον 10 τα set αυτά↩︎
Θυμηθείτε να με ρωτήσετε τι σημαίνουν οι υπόλοιπες μεταβλητές↩︎
Διαφορετικά δεν θα γίνει δεκτή και δεν θα βαθμολογηθείτε↩︎
Πάντα σύμφωνα με το είδος σας↩︎
Παραθέστε ΟΛΑ τα στοιχεία από τον σχετικό πίνακα↩︎
Delivered to you by Dr. Kostas Kougioumoutzis
kkougiou@aua.gr